clc;clear all;close all;
%main file to solve for the competitive equilibrium before the reform
global alpha betta depr sigmaC sigmaL gassetrat tauL tauK taumax lam
global phi popweight N govexp govery
global e_real omega kg

alpha = 0.394;   %capital share
betta = 0.96;   %discount factor
depr = 0.08;    %depreciation rate
sigmaC = 2;   %coefficient of relative risk aversion
sigmaL = 3;   %inverse Frisch elasticity

% tauL = 0.23;  %tax rates under the status quo policy 
% tauK = 0.57;
tauL=0.214 %Average of U.S. labour income tax rates taken from Trabandt and Uhlig 2012, period 2001-2010.
tauK=0.401  %Average of U.S. capital income tax with depreciation allowance taken from Trabandt and Uhlig 2012, period 2001-2010
taumax = tauK;  %max. capital tax rate
gassetrat = -0.668;  %initial government debt as a share of output

e_real=0.245

%HETEROGENOUS ECONOMY BEFORE THE REFORM

N=2

popweight = [.5 .5];  %population shares
phi = [1 1.1]'; %labor productivity
phi = phi/(popweight*phi); %normalize average labor productivity to 1
N = length(popweight);
lam=1/0.54; %c_j/c_1

options = optimset('MaxFunEvals',1000000,'MaxIter',1000,'TolFun',1e-8,'TolX',1e-8,'Disp','iter');
x0=[0.1  -0.2 2] %g kg k
lb=[0.01  -1 1.2]
ub=[0.4  0 2.5]
govery=0.2
[GE,info] = lsqnonlin('find_gkNoL',x0,lb,ub,options);
govexp = GE(1)
kg=GE(2)
kss=GE(3)
e_real=0.245

depr=alpha*e_real^(1-alpha)*kss^(alpha-1)-(1/betta-1)/(1-tauK); %the Euler at ss gives depr
css=(kss^alpha*e_real^(1-alpha)-govexp-depr*kss)/(popweight(1)+sum(popweight(2:N).*lam)); %the resource constraint gives c1
wss=(1-alpha)*kss^alpha*e_real^(-alpha);
%find l1 from optimality conditions of ratio of hours of different groups
%and using that average hours must match the data
l1ss=e_real;
k1ss=css^(-sigmaC)*(css-(1-tauL)*wss*phi(1)*l1ss)/(1-betta)/css^(-sigmaC)/(1+(alpha*(e_real/kss)^(1-alpha)-depr)*(1-tauK)); %lifetime budget constraint of group 1
k2ss=css^(-sigmaC)*(lam*css-(1-tauL)*wss*phi(2)*l1ss)/(1-betta)/css^(-sigmaC)/(1+(alpha*(e_real/kss)^(1-alpha)-depr)*(1-tauK));

cistst = [css lam*css]; %consumption of groups in the long run
listst = [l1ss l1ss]; %hours worked of groups in the long run
%utilities at the status quo
if (sigmaC==1)
    VsqNoL = 1/(1-betta)*log(cistst);
else
    VsqNoL = 1/(1-betta)*(cistst.^(1-sigmaC)-1)/(1-sigmaC);
end
VsqNoL

cistst_bm=cistst;

save('benchN2NoL.mat')

rss = (betta^(-1)-1)/(1-tauK)+depr; %r from Euler at the steady state
(rss/alpha)^(1/(1-alpha))*kss % express real total effective labor from r=F_k
govexp/kss^alpha/e_real^(1-alpha)
kss/kss^alpha/e_real^(1-alpha)

